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ABSTRACT 

The  objective  of  this  work  was  to  analyze  a  computer 
program  using  three  dimensional  quadratic  isoparametric 
finite  elements  for  structural  analysis.   Three  problems 
with  classical  solutions  were  run  with  various  mesh  sizes 
using  the  computer  program  being  tested.   The  data  computed 
was  then  extensively  analyzed,  and  compared  with  the 
classical  solutions.   The  analysis  of  a  fourth  problem  was 
continued  and  compared  with  results  obtained  in  an  earlier 
project . 
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I.   INTRODUCTION 

Before  the  advent  of  the  finite  element  technique  and 
the  digital  computer,  the  solutions  to  most  non-trivial 
elasto-static  problems  were  unobtainable.   Since  that  time 
many  finite  elements  have  been  devised  with  corresponding 
computer  programs  for  their  use. 

One  of  the  most  versatile  family  of  finite  elements  was 
developed  by  Professor  Zienkiewicz  and  his  co-workers  at 
the  University  of  Wales,  Swansea,  U.K.   This  family  of 
elements,  called  "ISOPARAMETRIC",   has  been  investigated  at 
the  Naval  Postgraduate  School  by  Professor  G.  Cantin,  and  a 
stress  analysis  program  called  TRISOP  was  written  using  a 
20  nodal  point  quadratic  element. 

Four  problems  with  known  solutions  will  be  solved  using 
TRISOP  and  analyzed  in  this  thesis.   From  this  analysis 
conclusions  will  be  made  as  to  the  effectiveness  and  accuracy 
of  TRISOP. 


II.   QUADRATIC  ISOPARAMETERIC  FINITE  ELEMENTS 

The  three-dimensional  element  discussed  in  this  chapter 
is  the  type  of  element  used  in  the  program  analyzed  by  this 
author.   No  attempt  will  be  made  to  fully  describe  this 
type  of  element,  but  just  to  give  some  insight  as  to  how 
the  element  is  constructed  and  designed. 

The  three-dimensional  quadratic  isoparameteric  finite 
element j  hereafter  referred  to  as  the  "element",  is  a  cube 
with  all  sides  two  units  in  length  in  it  undeformed,  non- 
dimensional  state.   The  element  has  20  nodal  points  and  is 
oriented  with  its  geometric  center  at  the  origin  of  its 
non-dimensionalized  coordinates,  (£,n,c),  and  transferred 
to  global,  (X,Y,Z),  coordinates  as  shown  in  Figure  1. 

A.   TRANSFORMATIONS 

Transformations  from  global  coordinates  to  non- 
dimensional  coordinates  are  accomplished  using  shape  functions 
as  shown  in  equations  1 

x(5,n,c)  =  NjL(?,n»c)x1 

yU,n,c)   =  N1(C,n,c)yjL  (D 

z(c,n,c)   =  Ni(c,n,c)zi 

where  N.(£,n,c)  are  shape  functions,  and  x.«,y.  >^^   are 
coordinates  of  the  nodal  points.   Similarly,  displacement 
transformations  are  made  using  equations  2, 


Figure  1.   Element  Transformation 
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u(e,n,0   =  N1(c,n,c)ui 
v(c,n,c)   =  Ni(c,n,c)v1 


w(£,n,c)   =   N1(c,n,c)wi 


(2) 


where  u,  v,  and  w  are  displacements  in  the  global  (X,Y,Z) 
reference  frame. 

The  transformation  of  a  volumetric  increment  from 
local  (c,n,c)  to  the  global  (X,Y,Z)  system  of  reference  is 
shown  in  equation  3.   The  Jacobian,  [J]  used  in  the 


dxdydz      =      det[J]   d^dndc 


(3) 


transformation  is  shown  in  equation  4. 


[J]   = 
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B.   SHAPE  FUNCTIONS 

For  an  element  with  nodes  numbered  as  in  Figure  2,  the 
shape  functions  obtained  from  Reference  1  are  shown  in 
equations  5. 

Corner  Nodes:   1,  3,  5,  7,  13,  15,  17,  and  19 
N±  =  (1/8)(1  +  £0)(1  +  n0)U  +  ^0)U0  +  %  +  ^0  -2)  (5) 

where  £   =  £ . £  and  £.  =  ±1;  similarly  for  n   and  c  . 
sO     1         1  o       o 
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Figure  2.   Element  Nodal  Point  Numbering  System 
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Midside  Nodes:   2,  6,  14,  and  18 
N±  =  (1/4X1  -  C2)(l  +  n0)(l  +  c0) 

Midside  Nodes:   4,  8,  16,  and  20  (5) 

Ni  =  (1/4) (1  -  n2)Cl  +  50)(1  +  ?Q) 

Midside  Nodes:   9,  10,  11,  and  12 
N±  =  (1/4  Xl  -  c2)(l  +  eo)(l  +  nQ) 

One  should  take  note  that  if  any  other  nodal  numbering 
system  with  respect  to  £  ,  n  ,  C  is  used,  the  nodal  point 
numbers  that  go  with  the  above  shape  functions  must  be 
adjusted  accordingly.   It  should  also  be  noted  that  the 
program  under  analysis  uses  the  system  mentioned  above. 

C.   COMPUTING  STRESS  AND  STRAIN 

Nodal  point  strains  are  computed  using  the  following 
relationship . 
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where  u,  v,  and  w  are  obtained  from  equation  2.  The  deri- 
vatives of  the  displacements  in  global  X,  Y,  Z  coordinates 
are  given  by: 
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where  [J]  is  the  Jacobian  matrix  (equation  4).   The  stresses 
are  computed  from  strains  by: 
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in  which:   G  = 


vE 


2(l+v)   '   A    (l+v)(l-2v) 


(9) 


E  is  Young's  modulus,  and  v  is  Poisson's  ratio 
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III.   TRISOP:   A  QUADRATIC  ISOPARAMETRIC 

FINITE  ELEMENT  PROGRAM 

The  stated  objective  of  this  work  is  to  evaluate  TRISOP, 
and  give  some  insight  into  its  practical  use.   Reference  1 
contains  the  theory  and  methods  used  in  the  program,  and 
Figure  3  is  a  simplified  flow  diagram. 

A.   INPUT  DATA 

The  input  data  required  by  TRISOP  are : 

1.  The  number  of  elements 

2.  The  total  number  of  nodal  points 

3.  The  number  of  different  materials 

4.  The  block  size  for  the  large  capacity  solver 

5.  The  number  of  nodal  points  with  boundary  conditions 

6.  The  number  of  nodal  points  with  concentrated  loads 

7.  Element  connectivity  for  each  element 

8.  The  coordinates  of  each  nodal  point 

9.  Young's  Modulus,  Poisson's  Ratio  for  each  material 

10.  Concentrated  loads 

11.  Boundary  conditions 

The  element  connectivity  is  a  correlation  between  the  over- 
all mesh  numbering  system  and  the  standard  numbering  system 
used  for  each  element. 

When  large  problems  are  to  be  solved,  preparing  the 
connectivity  and  coordinate  input  is  not  a  trivial  task. 
Also,  if  there  is  a  pressure  or  gravity  load,  the  computa- 
tions needed  to  convert  into  corresponding  nodal  point 
loads  can  be  extremely  time  consuming.   There  is,  however, 
a  mesh  generator  program  [2]  that  will  solve  this  problem. 
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With  a  minimum  input,  the  mesh  generator  will  compute  the 
element  connectivity,  nodal  point  coordinates  with  input  in 
cylindrical  or  rectangular  coordinates,  and  nodal  point 
loads  for  pressure  or  gravity  loads.   The  mesh  generator 
will  print  the  output,  draw  a  two  dimensional  picture  of 
the  mesh  and  punch  data  cards  for  TRISOP.   The  only  other 
input  to  TRISOP  that  can  require  more  than  a  few  cards  is 
the  boundary  conditions.   The  mesh  generator  will  not 
produce  boundary  condition  cards,  but  they  are  easily  and 
quickly  produced  by  hand. 

B.   TECHNIQUES  FOR  MOST  EFFICIENT  USE  OF  TRISOP 

TRISOP,  in  its  present  form,  has  a  constant  core  storage 
requirement.   The  two  variables  with  the  size  of  the  problem 
are  disk  storage  requirements  and  running  time.   To  make  the 
most  efficient  use  of  TRISOP,  it  is  essential  to  take  advan- 
tage of  symmetry  whenever  possible.   It  is  also  essential 
to  design  the  problem  mesh  in  order  to  reduce  the  half-band 
width  of  the  resulting  system  of  equations  to  a  minimum. 
The  half-band  width  is  a  function  of  the  difference 
between  the  highest  and  lowest  nodal  point  number  in  any 
element.   To  obtain  the  smallest  band  width  one  must  start 
numbering   on  the  face  having  the  least  number  of  elements 
along  the  side  having  the  smallest  number  of  elements  as 
shown  in  Figure  4 . 

TRISOP  does  not  give  exact  answers  to  a  problem,  but 
will  converge  asymptotically  and  monotonically  to  an  exact 
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solution  with  uniform  mesh  refinement.   To  insure  reliable 
results  to  a  problem,  a  convergence  study  must  be  made. 

The  convergence  technique  used  here  consists  of  plotting 

p 

the  displacement  of  a  convenient  nodal  point,  versus  1/N 

where  N  is  the  number  of  elements  in  the  mesh.   If  three 
points  plot  in  a  straight  line,  the  extrapolation  to  the 
origin  is  justified.   This  convergence  technique  is  an 
adaptation  of  a  technique  developed  by  L.  P.  Richardson  [3] 

Richardson's  technique  was  developed  for  extrapolating 
the  results  of  central  finite  difference  approximations 
where  the  truncation  error  is  of  the  order  (h  ),  and  h  is 
the  finite  difference  interval.   When  the  truncation  error 
is  of  the  order  (h  )  the  extrapolated  value  is  found  by: 


"extrap 


2       2 
Xph   -  x.hp 

hl   "  h2 


(10) 


This  formula  is  plotted  in  Figure  5. 


extrap 


h  2  u  2     .2 
h3   h2     hl 


Figure  5.   Richardson's  Technique 
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If  a  third  value,  x~,  is  found  with  a  third  interval,  h  , 
and  x^  plots  on  the  extrapolation  line,  the  extrapolated 
value  is  presumed  to  be  exact  [3]. 

There  is  no  guarantee  that  the  error  in  TRISOP  is  of 

2 
order  1/N  ,  but  this  method  of  ascertaining  convergence 

was  adopted  after  some  numerical  experimentation.   Since 
it  is  assumed  that  the  solution  approaches  the  exact  solu- 
tion asymptotically,  if  three  deflections  plot  in  a  straight 
line  it  is  assumed  that  extrapolation  is  valid.   If  they 
don't  plot  in  a  straight  line  there  is  no  way  of  predicting 
the  accuracy  of  the  extrapolated  result.   This  technique 
has  worked  quite  well  with  many  problems,  but  a  better 
means  of  determining  convergence  is  needed. 

C.   IMPROVEMENTS  MADE  TO  TRISOP 

Gaussian  integration  is  used  extensively  in  the  compu- 
tation of  the  element  stiffness.   Initially,  four  Gauss 
points  in  each  of  the  three  directions  of  a  coordinate 
system  were  used  in  the  solution.   The  subroutine  in  which 
this  was  accomplished  was  called  CUE4 .   During  a  visit  by 
Professor  0.  C.  Zienkiewicz  this  author  was  told  that  using 
two  Gauss  points  in  the  integration  improved  the  solution. 
Subroutine  CUB2,  a  two  point  integration  subroutine,  was 
substituted  for  CUB4  in  TRISOP.   This  change  yielded  better 
results  with  a  coarser  mesh  for  all  problems  tested.   The 
change  from  CUE4  to  CUB2  also  reduces  significantly  the 
integration  CPU  time  for  the  calculation  of  element 
stiffness. 


20 


If  an  element  is  deformed  into  an  extreme  shape  such 
as  the  element  in  Figure  6,  the  Jacobian  becomes  singular 
at  points  similar  to  3,  10,  and  15.   Since  equation  7 
which  uses  the  inverse  of  the  Jacobian  is  used  in  computing 
stress  and  strain  values,  the  computer  algorithm  fails  at 
points  3,  10,  and  15.   To  eliminate  this  singularity, 
TRISOP  was  modified  to  displace  each  node  by  a  small  dis- 
tance away  from  its  actual  location  at  the  time  the 
Jacobian  is  formed  for  that  node.   The  flow  of  computations 
is  then  uninterrupted,  and  the  results  for  all  other  nodes 
are  unaffected.   The  stresses  and  strains  computed  for  nodes 
with  a  singularity  are  meaningless  and  should  be  eliminated 
from  any  further  considerations. 
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Figure  6.   Degenerate  Element 
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IV.   SOLUTIONS  TO  CLASSICAL  PROBLEMS  USING  TRISOP 

Four  classical  problems  were  considered  in  the  analysis 

of  TRISOP.   The  four  problems  were  a  simply  supported  beam, 

a  pinched  disk,  the  Boussinesq  problem,  and  a  pinched 
cylinder. 

A.   SIMPLY  SUPPORTED  BEAM 

The  simply  supported  beam  shown  in  Figure  7  was  analyzed 

1.   Classical  Solution 

The  classical  solution  used  was  an  Airy  stress 
function  developed  using  elasticity  theory  [4].   The  Airy 
stress  function  for  the  problem  under  consideration  is 
shown  in  equation  11. 


*  -  -  l  x2  ♦  §a  x2y  -  -aj  xV  ♦  &<4  -  f)  73  +  -*?  y5     (id 

*  oc       8c  J        oc  c     °  40cJ 


Equation  11  satisfies  equation  12  as  is  required 


i-4  +  2    p  ^2  +  ^"T   =   °  (12) 

ax     ax  ay    ay 


for  a  valid  Airy  stress  function  in  the  absence  of  body 

forces.   Stress  values  are  found  from  equation  11  using 
equation  13. 

2  2  2 

°*~-7?     '      °'~I?     '       x*>"w  (13) 
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Evaluating  these  relations  yields  equation  lk 


a,  -  £  (L2-x2)y  +  &  (  f  y3  -  §  o*y> 


ay  =  -  A  (  i  y3  -  °2y  +  f  °3)  t1^ 


2    2 
xy     21 


s  -  3t  (c  -  y  )x 


2c3 
where  I  =  — ^—  is  the  moment  of  inertia  of  the  beam.   Since 

there  are  no  body  forces,  such  as  gravity  loading,  this  is 
the  solution  to  both  the  plane  stress  and  plane  strain 
problems.   The  stress  component,  sigmaxis  an  exact  solu- 
tion to  the  problem  only  if  the  distributed  normal  force 
shown  in  equation  15  is  applied 

*   -    ±  i  ■%  (  I  y3  -  I  c2y)  (i5) 

c 

to  each  end  of  the  beam.   Since  this  force  distribution 
has  zero  resultant  force,  and  zero  resultant  moment,  it  can 
be  concluded  by  Saint-Venant ' s  Principle  [4]  that  sigma x 
is  exact  at  some  distance  from  the  ends  of  the  beam. 
Saint -Venant ' s  Principle  assumes  that  localized  forces  in 
static  equilibrium  will  give  rise  to  localized  stresses 
and  strains. 

The  displacement  of  the  center  of  the  beam  is 
found  using  equation  16.   This  displacement  is  greater 

«■   ^tr^1  +  T  ^4  +  f)]  (16) 
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than  the  displacement  found  using  elementary  theory  because 
there  the  assumption  is  made  that  cross  sections  of  the  beam 
remain  plane  during  bending. 
2.   TRISOP  Formulation 

Considering  the  problem  in  Figure  7,  there  is  a 
plane  of  symmetry  parallel  to  the  XY  plane  passing  through 
the  Z  axis  at  .5  inches.   All  the  points  in  this  plane  of 
symmetry  were  given  a  geometrical  constraint  which  set  the 
displacement  component  in  the  Z  direction  equal  to  zero. 
Only  half  of  the  beam  from  0  <_  Z  <_  .  5  was  considered.   For 
the  nodal  points  in  the  XZ  plane  at  X  =  ±L,  the  displace- 
ment component  in  the  Y  direction  was  set  equal  to  zero. 
Finally,  for  the  nodal  points  on  the  YZ  plane  the  displace- 
ment component  in  the  X  direction  was  set  equal  to  zero. 
The  distributed  load  of  1000  psi  is  modeled  using  a  consist- 
ent load  vector  [5].   This  vector  dictates  that  the  force 
on  the  face  of  each  loaded  element  be  divided  as  shown  in 
Figure  8. 


F  =  qA 


Figure  8.  Consistent  Load  Vector 
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3 .   Results 

A  convergence  study  of  the  deflections  at  the 
center  of  the  beam  is  shown  in  Figure  9  using  the  meshes 
shown  in  Figure  10.   The  results  indicate  that  extrapola- 
tion to  a  mesh  composed  of  an  infinity  of  elements  is 
justified.   The  classical  value  for  the  deflection  at  the 
center  of  the  beam  is  0.0180  inches.   The  extrapolated 
value  given  by  TRISOP  was  0.0188  inches ,  a  value  .0008 
inches  greater  than  the  classical  value.   The  classical 
solution,  being  a  two  dimensional  solution,  assumes  that 
there  is  no  deflection  perpendicular  to  the  XY  plane. 
This  makes  the  classical  solution  stiffer  than  the  three 
dimensional  solution.   Figure  11  shows  the  displacement 
component  in  the  Z  direction  in  the  XY  plane  given  by  TRISOP 
When  TRISOP  results  for  sigma  x,  sigma  y,  and  tau  xy 
were  compared  with  the  classical  solution  using  nodes  on 
the  mesh  face  the  values  given  by  sigma x  were  within  at 
least  .1  percent  of  the  classical  value  for  all  nodal 
points.   On  the  other  hand,  the  values  for  sigma  y  and 
tau  xy  were  not  nearly  as  good.   It  was  noted,  however,  that 
the  best  results  were  obtained  from  interior  nodes  in  the 
mesh.   A  8x3x2  mesh  was  then  analyzed  using  data  from  the 
mid  plane  parallel  to  the  XY  plane.   The  results  from  this 
analysis  are  shown  in  Figures  12,  13,  and  1 4 .   Since  the 
beam  is  symmetrical  about  the  origin  along  the  X  axis, 
only  the  left  end  of  the  beam,  -L  <_  X  <_   0 ,  is  shown. 


27 


CM 


to 

<D 

E 


X 

CM 

X 

00 


^3- 

4    O 

O 


JC 

to 

CD 


X 
CO 
X 

00 


to 


X 

«3- 

X 

00 


4-> 

c 

0) 

E 
CD 

o 

(O 

*D- 
to 

•r- 

TO 

h- 

Ift 

CM 

o 


to 
o 

to  CD 

to  3 

rO  i— 

i—  fO 

a  > 


CO 

o 
o 


CM 

o 
o 


o 
o 


o 

CM 


in 
o 


>i 

-J 

-p 
co 

CD 

o 

0) 

b£) 

u 

> 
c 
o 
o 

s 

<u 
P3 

-p 
o 

co 

>s 

rH 

6 
•H 
CO 


ON 

<D 
U 

w 

•H 


o 
o 


28 


T r 


8x2x1  mesh 


8x3x1   mesh 


8x4x1  mesh 


y^        /        y^        y^       y^      y^        y^        y^ 

y 

> 

y^ 
y 

! 

i 

• 

J  __ 

yS            // 

^S»                  */* 

y^ 

// 

y 

1 

A 

A 

i 

i 

j 

A 

i 
1 

- 

i 

Figure    10.      Simply   Supported   Beam  Meshes 
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Displacements  x  10  J   on  XY  Plane  of  the 
Simply  Supported  Beam ■ 

Figure  11 
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Figure  12.   Simply  Supported  Beam 
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Figure  13.   Simply  Supported  Beam 
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X=-L 


Classical  Solution 


X=0 


0 

The  lower  half  is  the  image  of  the  upper 
half.  The  right  end  is  the  image  of  the 
left  end  with  reversed  signs. 
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Figure  14.   Simply  Supported  Beam 


rxy  x  10  psi 


33 


4.   Conclusions 

The  results  of  this  analysis  indicate  that  much 
more  reliable  data  is  obtained  from  interior  nodes.   The 
poor  results  obtained  for  sigma  y  and  tau  xy  on  the  mesh  faces 
could  indicate  that  the  values  obtained  under  a  consistent 
load  vector  are  better  in  planes  where  the  loaded  nodal 
points  have  four  common  elements  as  in  the  mid  plane  of  the 
8x3x2  mesh. 

B.   PINCHED  DISK 

A  pinched  disk  with  two  equal  and  opposite  forces 
acting  on  the  diameter  as  shown  in  Figure  15  was  analyzed. 

1.   Classical  Solution 

A  classical  two  dimensional  solution  was  devised 
using  an  Airy  stress  function  by  H.  Hertz,  and  is  discussed 
in  detail  in  Reference  4.   The  Airy  stress  function, 
equation  17,  gives  the  stress  in  terms  of  e  and  r  for  each 
load  as  shown  in  Figure  16.   The  radial  stress,  the  only 

<t>  =  -  re  sin  e  (17) 

non  zero  stress,  is  shown  in  Equation  18.   Equation  18 

1    3<(>         1     92<}>  2P   cos    e  /noN 

°r  =  ?^    7  J?  "  ~—  (18) 

leads  to  the  presence  of  an  isotropic  state  of  compression 
of  intensity  2P/-rrd  all  around  the  radial  surface  of  the 
disk.   To  free  the  boundary  of  this  unwanted  stress  an 
isotropic  tension  equal  to  2P/ird  is  added  to  the  disk. 
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400  K 


Thickness  =  1" 
E  =  30  x  TO6  psi 
v  =  0.3 


Figure    15.      Pinched   Disk 
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Figure  16.   Pinched  Disk  Radial  Stress 
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After  converting  sigma  r-,  and  sigma  r^  into  X,Y  coordinates 

and  adding  the  stresses  from  each  load  together,  equations 
19  result. 

2  2      2P 

a    =   -a-  sin  6,  -  a   sin   6~  +  — r 
x       r,  r? 

2  2       2P 

a    =   -a   cos  6,-o   cos   e0  +  — r         (19) 
y       r-,  2       2    ird 


t     =   -a   sin  e,  cos  B_  +  o   sin  90  cos  e0 
xy       r1     1      1    r2     2      2 


2.  TRISOP  Formulation 

The  problem  was  formulated  for  TRISOP  as  shown  in 
Figure  17  using  half  the  thickness  of  the  disk.   The 
boundary  conditions  were  w  =  0  in  the  XY  plane,  u  =  0  in 
the  YZ  plane,  and  v  =  0  in  the  XZ  plane.   Figure  18  shows 
one  face  of  the  8x8x1  mesh.   To  model  the  load,  a  consistent 
load  vector  was  again  used  [5]  with  values  as  shown  in 
Figure  17. 

3.  Results 

The  deflections  on  the  radial  surface  of  the  disk 
at  a  distance  of  0.975  inches  from  the  load  were  0.01130 
inches  for  the  4x4x1  mesh,  and  0.01133  inches  for  the 
8x8x1  mesh.   Due  to  the  small  difference  in  these  two 
deflections,  it  was  assumed  that  the  solution  could  be 
extrapolated  with  no  significant  error.   The  convergence 
plot  is  shown  in  Figure  19. 
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=  100  K 


w  =  0 


u  =  0 
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Figure  17.   Pinched  Disk,  TRISOP  Formulation 
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Figure  18.   Pinched  Disk  8x8x1  Mesh 
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Figure  19.   Pinched  Disk  Convergence  Study 
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Contour  graphs  were  drawn  using  data  from  the  8x8x1 
mesh  to  compare  TRISOP  data  with  the  classical  solution 
obtained  by  Carlos  Felippa  [6].   Two  graphs  were  drawn,  one 
using  mid  side  node  data,  and  one  using  corner  node  data 
with  the  results  shown  in  Figure  20.   The  plots  for  sigma  y 
and  tau  xy  using  corner  and  mid  side  nodes  are  shown  in 
Figures  21  and  22.   All  TRISOP  plots  shown  used  data  from 
the  constrained  (XY)  plane.   However,  the  data  on  the  free 
plane  (Z  =  .5  inches),  and  the  mid  plane  (Z  =  .25  inches) 
was  virtually  the  same  as  that  in  the  corresponding  nodal 
points  in  the  XY  plane. 

k .   Conclusions 

The  TRISOP  solution  agrees  with  the  classical  solu- 
tion for  sigma  y,  tau  xy,  and  the  mid  side  node  data  for 
sigma  x.   The  results  for  sigma  x  using  corner  node  data 
indicates  a  much  higher  state  of  stress  in  the  X  direction 
directly  under  the  load.   It  should  be  noted,  however,  that 
there  are  only  two  nodes  that  have  values  that  are  greatly 
in  error,  and  one  of  those  is  directly  under  the  load. 

C.   BOUSSINESQ  PROBLEM 

The  Boussinesq  problem  consists  of  one  concentrated 
load  normal  to  the  surface  of  a  semi-infinite  solid  as 
shown  in  Figure  23. 
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TRISOP  Solution 

using  corner  node  data 


25  20  15  10  5 


TRISOP  Solution 

using  mid  side  node  data 


Classical  Solution 


Figure  20.   Pinched  Disk  a      x  103  psi 
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TRISOP  Solution 


Classical  Solution 


Figure  21.   Pinched  Disk  a      x  10   psi 
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TRISOP  Solution 


Classical  Solution 


Figure  22.   Pinched  Disk   x    x  10   psi 
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Figure  23.   Boussinesq  Problem 
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1.   Classical  Solution 

The  classical  solution  to  the  problem  of  a 
concentrated  load  normal  to  the  surface  of  a  semi-infinite 
solid  was  solved  by  J.  Boussinesq  [4]  and  yields  the  results 
shown  in  equations  20.   The  displacement  in  the  Z  direction 
is  given  in  equation  21. 


c 


a 


»  =   oT  i(l-2v)[— o  ~  ~r(r  +z  )     J~3r  z(r  +z  )     } 

r    r 


3P   3/  2^  2,-5/2 


z       2tt 


(20) 


P  ,,  «  v,  1     z  ,  2^  2,-1/2     ,  2^  2x-3/2, 
r    r 

3P    2,  2^  2,-5/2 
Trz  =   "  27  rz  (r  +z  } 


w    =   «£=■  [(l+v)z2(r2+z2)~3/2  +  2 (1-v2 ) (r2+z2)"1/2]     (21) 
d.~nhj 


In  equations  20  and  21,  E  is  Young's  Modulus,  and  v  is 
Poisson's  ratio.   The  stresses  are  in  cylindrical  coordi- 
nates with  the  axis  of  symmetry  being  the  line  of  action 
of  the  load. 

2.   TRISOP  Formulation 

The  TRISOP  formulation  took  advantage  of  the  double 
symmetry  of  the  problem,  and  is  shown  in  Figure  24.   This 
figure  also  shows  the  dimensions  of  the  3x3x3  mesh.   The 
2x2x2,  and  4x4x4  meshes  were  similarly  constructed,  and 
Figure  25  shows  one  face  of  these  meshes  to  indicate  the 
element  dimensions. 
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Figure  2k.      Boussinesq  Problem  TRISOP  Formulation 
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Figure  25.   Boussinesq  Meshes 
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3.  Results 

As  evidenced  by  equation  21,  the  deflection  at  the 
load  is  not  bounded.   After  examining  deflections  at  other 
points  in  the  three  meshes,  it  could  not  be  determined  if 
extrapolation  to  a  converged  solution  would  be  valid  for 
the  meshes  used.   Figure  26  shows  the  solution  of  the 
deflections  in  the  Z  direction  from  r  =  0  to  r  =  100,  at 
distances  of  5,  10,  and  20  inches  from  the  surface  using 
nodal  point  deflections  from  the  kxkxk   mesh.   The  TRISOP 
deflections  follow  the  general  contour  of  the  classical 
solution,  but  are  not  as  large. 

Figures  27,  28,  29,  and  30  show  contour  graphs  of 
sigma  r,.  sigma  e,  sigma z  and  tau  rz  using  classical  and  TRISOP 
data.   A  block  twenty  inches  on  a  side  is  used  for  these 
graphs,  because  the  area  close  to  the  load  is  of  prime 
interest.   The  data  used  in  generating  the  TRISOP  graphs 
came  from  interior  nodal  points  in  the  mesh.   Since  TRISOP 
computes  results  in  rectangular  coordinates,  and  the 
Boussinesq  solution  is  in  cylindrical  coordinates,  it  was 
necessary  to  transform  stress  data  from  the  interior  nodes 
to  cylindrical  coordinates  for  meaningful  comparison  with 
the  classical  solution. 

4 .  Conclusions 

The  discrepancy  between  the  classical  solution  for 
displacements  and  the  TRISOP  solution  is  caused  by  the 
problem  formulation  used.   In  the  real  problem  the  deflec- 
tions at  Z  =  100  inches  are  not  zero.   If  the  boundary 
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Figure    26.      Boussinesq  Deflections 


50 


TRISOP 
Solution 


a      x  10  psi  with  r  and  z  from  0  to  20  inches 
z      F 


Classical 
Solution 


Figure  27.   Boussinesq  Problem 
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Figure  28-.   Boussinesq  Problem 
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Figure  29.   Boussinesq  Problem 


53 


TRISOP 
Solution 


T   x  10  psi  with  r  and  z  from  0  to  20  inches 
lrz      ^ 


Classical 
Solution 


Figure  30.   Boussinesq  Problem 
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conditions  at  the  base  of  the  mesh  could  have  modeled  the 
expected  deflections,  it  is  believed  that  the  deflections 
elsewhere  in  the  mesh  would  have  been  nearer  to  the  expected 
values.   The  contour  graphs  indicate  that  TRISOP  generates 
a  solution  in  close  agreement  with  the  classical  solution. 

D.   PINCHED  CYLINDER 

The  pinched  cylinder  consists  of  a  thin  shell  cylinder 
with  a  concentrated  radial  load  as  shown  in  Figure  31. 
D.  E.  Hanson  [2]  solved  this  problem  using  TRISOP  with  four 
Gauss  point  integration,  and  compared  his  results  to  a 
study  mady  be  G.  Cantin  [7]  using  thin  shell  elements.   The 
objective  of  this  study  was  to  determine  if  the  three 
dimensional  solution  would  give  comparable  results.   Hanson 
made  runs  with  meshes  up  to  a  1x10x10,  and  obtained 
discouraging  results.   After  discovering  that  two  Gauss 
point  integration  might  yield  better  results,  Hanson's 
1x2x2  and  1x^x4  meshes  were  rerun  using  two  point  integra- 
tion.  Table  1  shows  a  tabulation  of  Hanson's  and  Cantin's 
results  including  the  reruns  using  two  point  integration. 
Figure  32  is  a  convergence  study  of  the  two  meshes  run  with 
two  point  integration.   The  extrapolated  result  agrees  with 
the  fully  converged  solution  given  by  Cantin. 
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Figure    31.      Pinched   Cylinder 
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TABLE  I 
Displacement  of  a  Pinched  Cylinder  at  the  Applied  Load 


Cantin 


Hanson 


No. 

of 

Mesh 

Displac . , 

Eq. 

in. 

36 

1x2 

-0.0921 

48 

1x3 

-0.1072 

60 

1x4 

-0.1099 

54 

2x2 

-0.0931 

72 

2x4 

-0.1085 

90 

2x4 

-0.1113 

150 

4x4 

-0.1126 

294 

6x6 

-0.1137 

486 

8x8 

-0.1139 

726 

10x10 

-0.1139 

No. 

of 

Mesh 

Displac . , 

Eq. 

in . 

153 

1x2x2 

-0.00177 

465 

1x4x4 

-0.00375 

945 

1x6x6 

-0.03002 

1593 

1x8x8 

-0.08702 

2409 

1x10x10 

-0.10057 

Hanson  with  Two  Point  Integration 


No. 

of 

Mesh 

Displac 

Eq. 

in. 

153 

465 


1x2x2 
1x4x4 


-0.1043 
-0.1127 
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V.   CONCLUSIONS  AND  RECOMMENDATIONS 

The  analysis  of  the  problems  considered  indicates  that 
TRISOP  can  be  expected  to  give  accurate  results.   However 
it  was  found  in  the  simply  supported  beam  and  Boussinesq 
problems  that  data  obtained  for  external  surfaces  were 
inaccurate.   It  would  therefore  be  wise  to  avoid  using 
data  from  surface  nodes.   If  information  near  the  surface 
of  a  solid  is  desired,  a  possible  solution  would  be  to  have 
a  very  thin  plane  of  elements  at  the  surface,  and  use  the 
data  generated  by  the  elements  just  below  the  surface. 
Another  possible  solution  would  be  to  modify  the  program 
to  compute  stress  and  strain  at  points  other  than  nodal 
points  such  as  the  Gauss  integration  points  to  see  if 
better  results  are  obtained.   However,  this  would  increase 
the  complexity  of  the  program,  because  the  coordinates  of 
the  points  to  be  used  would  have  to  be  specified  or  calcu- 
lated and  exhibited  together  with  the  stress  values. 

Concurrent  with  this  research  LCDR  E.  Leonidas  was 
making  improvements  to  TRISOP  [8]  that  included  reducing 
the  program  run  time,  and  generalizing  the  boundary  condi- 
tions to  make  it  possible  to  specify  boundary  displacements. 
This  last  improvement  will  make  it  possible  to  model 
problems  such  as  the  Boussinesq  problem  more  accurately. 

TRISOP  has  the  capability  of  receiving  input  mesh  data 
in  cylindrical  coordinates,  but  this  data  is  converted,  and 
the  output  is  given  in  rectangular  coordinates.   For  problems 
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like  the  Boussinesq  problem,  output  data  in  cylindrical 
coordinates  would  be  a  great  help.   TRISOP  could  be  modified 
to  accomplish  this,  or  possibly  another  similar  program 
could  be  produced  to  solve  problems  in  cylindrical  and 
spherical  coordinates. 

The  change  from  four  point  to  two  point  integration 
greatly  reduces  the  computer  time  needed  to  solve  a  problem. 
This  change  also  seems  to  give  a  more  accurate  solution  as 
shown  by  the  results  of  the  pinched  cylinder  analysis.   The 
computer  time  needed  to  solve  large  problems,  however,  can 
be  excessive  and  is  a  very  real  limitation.   Table  II  shows 
the  CPU  time  for  various  problems  run  with  both  CUB2  and 
CUB4  on  an  IBM  360  computer. 

The  largest  shortcoming  of  TRISOP  is  that  the  excessive 
amount  of  data  produced  is  very  time  consuming  and  tedious 
to  analyze.   There  are  two  possible  solutions  to  this 
problem.   The  first  would  be  to  modify  the  program  so  that 
the  output  would  be  placed  on  contour  graphs  by  the  computer 
The  problem  with  this  solution  is  that  if  the  graphs  do  not 
plot,  or  graphs  other  than  those  asked  for  are  needed  after 
the  run,  the  program  would  have  to  be  rerun.   In  the  case 
of  a  large  problem,  the  cost  in  computer  time  could  be 
excessive.   What  appears  to  be  a  far  better  solution  would 
be  to  have  the  output  placed  on  magnetic  tape.   The  data 
could  then  be  thoroughly  analyzed  using  the  computer  with 
a  minimum  expenditure  of  computer  time. 
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TABLE  II 
COMPUTER  TIMES  FOR  TEST  PROBLEMS  RUN 


Simply  Supported  Beam 

1x2x8  CUB4 
1x3x8  CUB4 
1x4x8  CUB4 
2x3x8  CUB4 
1x2x8  CUB2 
1x3x8  CUB2 
1x4x8  CUB2 
2x3x8  OUB2 

Pinched  Disk 

1x4x4  CUB4 
1x5x5  CUB4 
1x6x6  CUB4 
1x8x8  CUB4 
1x4x4  CUB2 
1x5x5  CUB2 
1x6x6  CUB2 
1x8x8  CUB2 

Boussinesq  Problem 

2x2x2  CUB2 
3x3x3  CUB2 
4x4x4  CUB2 

Pinched  Cylinder 

1x4x4  CUB4 
1x4x4  CUB2 


Overall 

PSTF 

Merge 

Solve 

Time 

(sec) 

(sec) 

(sec) 

(min,  sec) 

310.38 

120.77 

381.27 

13 

53.58 

462.31 

174.73 

539.26 

20 

06.58 

619.38 

353.28 

1116.13 

35 

23.07 

897.62 

507.29 

1357.63 

47 

44.98 

48.49 

118.35 

386.38 

9 

33.73 

72.01 

174.58 

534.60 

13 

31.32 

94.47 

347.45 

1112.54 

26 

29.13 

135.66 

507.26 

1357.67 

35 

01.59 

315.05 

175.82 

522.73 

17 

09.66 

470.96 

285.88 

835.18 

26 

57.49 

675.46 

566.05 

1694.84 

49 

25.42 

1195.18 

1433.38 

4217.03 

115 

33.68 

47.38 

173.84 

502.34 

12 

19.66 

72.47 

288.26 

832.60 

20 

18.46 

102.60 

566.32 

1688.59 

39 

44.81 

179.07 

1433.12 

4216.28 

98 

34.08 

25.19 

90.43 

249.43 

6 

12.41 

76.61 

432.22 

1062.09 

26 

29.91 

178.31 

1853.68 

4309.12 

107 

07-95 

317.60 

152.09 

646.10 

18 

53.85 

44.16 

147.75 

497.40 

12 

42.00 
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